#include <fstream>
#include <sstream>
#include <cmath>
#include <boost/format.hpp>
#include "CPPDiffusion2DADI.h"

void PrepareInitialValues(double *, CPPDiffusion2DADI *);
void Print(CPPDiffusion2DADI *, double *, const char *);

int main()
{
	int endCount = 100;

	CPPDiffusion2DADI * sim = new CPPDiffusion2DADI();
	sim->BoundaryCondition(0.0);
	sim->NumberOfNodesInX(101);
	sim->NumberOfNodesInY(101);
	sim->DeltaX(1.0 / 100);
	sim->DeltaY(1.0 / 100);
	sim->DeltaT(0.01);
	sim->Alpha(0.005);

	double * u = new double[101 * 101];
	PrepareInitialValues(u, sim);

	for (int i = 0; i <= endCount; i++) {
		double timeInMiliseconds = i * sim->DeltaT() * 1000.0;
		std::stringstream buf;
		buf << boost::format("c:/nsworkspace/cppdiffusion2dadi_%05d.txt") 
			% timeInMiliseconds;
		Print(sim, u, buf.str().c_str()); 

		sim->NextStep(u);
	}
	return 0;
}

void PrepareInitialValues(double * u, CPPDiffusion2DADI * sim) {
	for (int j = 0; j < sim->NumberOfNodesInY(); j++) {
		for (int i = 0; i < sim->NumberOfNodesInX(); i++) {
			double x = i * sim->DeltaX();
			double y = j * sim->DeltaY();

			u[j * sim->NumberOfNodesInX() + i] = 
				std::exp(-0.5 * std::pow((x - 0.5) / 0.1, 2.0) 
				         -0.5 * std::pow((y - 0.5) / 0.1, 2.0));

			if (i == 0 || i == sim->NumberOfNodesInX() - 1 ||
				j == 0 || j == sim->NumberOfNodesInY() - 1) {
					u[j * sim->NumberOfNodesInX() + i] = 0.0;
			}
		}
	}
}

void Print(CPPDiffusion2DADI * sim, double * u, const char * outfile) {
	std::ofstream writer;		
	writer.open(outfile, std::ios::out);

	for (int i = 0; i < sim->NumberOfNodesInX(); i++) {
		for (int j = 0; j < sim->NumberOfNodesInY(); j++) {
			double x = i * sim->DeltaX();
			double y = j * sim->DeltaY();
			writer << x << " " << y << " ";
			writer << u[j * sim->NumberOfNodesInX() + i] << std::endl;
		}
		writer << std::endl;;
	}
	writer.close();
}


